package algorithm.sphericharmonics;

import algorithm.fft.Complex;

public class SphericalHarmonics extends MathFunct{
	static final float TWOPI=2*PI; 
	public static ComplexSphere y0(int m,float theta, float phi){
		if (m==0){
			return new ComplexSphere(sqrt(1/PI)/2,0); 
		}
		throw new RuntimeException("illegal params"); 
	}
	public static ComplexSphere yfunct(int l,int m, float theta,float phi){
		float f=1;
		int mm=m; 
		if (m<0){
			if (m%2==0){
				f=1; 
			}else {
				f=-1;  
			}
			mm=-m;
			phi=-phi; 
		}
		ComplexSphere c=null; 
		switch(l){
		case 0:
			c=y0(mm,theta,phi);
			break; 
		case 1:
			c=y1(mm,theta,phi);
			break; 
		case 2:
			c=y2(mm,theta,phi);
			break; 
		case 3:
			c=y3(mm,theta,phi);
			break; 
		case 4:
			c=y4(mm,theta,phi);
			break; 
		}
		if (c!=null){
			c.radius=c.radius*f;
			if (c.radius<0){
				c.angle+=PI;
				c.radius=-c.radius; 
			}
			return c; 
		}
		throw new RuntimeException("not implemented");
	}
	public static ComplexSphere y1(int m, float theta, float phi){
		switch(m){
			case -1:
				return new ComplexSphere(sqrt(3/TWOPI)/2*sin(theta),-phi); 
			case 0:
				return new ComplexSphere(sqrt(3/PI)/2*cos(theta),0);
			case  1:
				return new ComplexSphere(-sqrt(3/TWOPI)/2*sin(theta),phi);
		}
		throw new RuntimeException("illegal params");
	}
	public static ComplexSphere y2(int m,float theta,float phi){
		switch(m){
		case -2:
			return new ComplexSphere(sqrt(15/TWOPI)/4*sin(theta)*sin(theta),-2*phi);
		case -1:
			return new ComplexSphere(sqrt(15/TWOPI)/2*sin(theta)*cos(theta),-phi);
		case 0:
			return new ComplexSphere(sqrt(5/PI)*(3*cos(theta)*cos(theta)-1)/4.0f,0);
		case 1:
			return new ComplexSphere(-sqrt(15/TWOPI)/2*sin(theta)*cos(theta),phi);
		case 2:
			return new ComplexSphere(sqrt(15/TWOPI)/4*sin(theta)*sin(theta),2*phi);
		}
		throw new RuntimeException("illegal params");
	}
	public static ComplexSphere y3(int m, float theta, float phi){
		switch(m){
		case 0:
			return new ComplexSphere(sqrt(7/PI)/4*(5*cos(theta)*cos(theta)*cos(theta)-3*cos(theta)),0); 
		case 1:
			return new ComplexSphere(-sqrt(21/PI)/8*sin(theta)*(5*cos(theta)*cos(theta)-1),phi); 
		case 2:
			return new ComplexSphere(sqrt(105/TWOPI)/4*sin(theta)*sin(theta)*cos(theta),2*phi); 
		case 3:
			return new ComplexSphere(-sqrt(35/PI)/8*sin(theta)*sin(theta)*sin(theta),3*phi); 
		}
		throw new RuntimeException("illegal params");
	}
	public static ComplexSphere y4(int m, float theta, float phi){
		switch(m){
		case 0: 
			return new ComplexSphere((3.0f/16)*sqrt(1/PI)*(35*pow(cos(theta),4)-30*cos(theta)*cos(theta)+3),0); 
		case 1:
			return new ComplexSphere((-3.0f/8)*sqrt(5/PI)
					*sin(theta)*(7*cos(theta)*cos(theta)*cos(theta)-3*cos(theta)),phi); 
		case 2:
			return new ComplexSphere((3.0f/8)*sqrt(5/TWOPI)*sin(theta)*sin(theta)*(7*cos(theta)*cos(theta)-1),2*phi);
		case 3:
			return new ComplexSphere((-3/8.0f)*sqrt(35/PI)*sin(theta)*sin(theta)*sin(theta)*cos(theta),3*phi);
		case 4:
			return new ComplexSphere((3/16.0f)*sqrt(35/TWOPI)*sin(theta)*sin(theta)*sin(theta)*sin(theta),4*phi);
		}
		throw new RuntimeException("illegal params");
	}
	public static void test(float theta,float phi){
		for (int l=0; l<=4; l++){
			for (int m=-l; m<=l; m++){
				ComplexSphere c1=SphericalHarmonics.yfunct(l, m, theta, phi); 
				ComplexSphere c2=ScaledAssociatedLegenedrePolynomialFunct.yfunction(l, m, theta, phi);
				//Complex c=c1.toComplex().minus(c2.toComplex());
				double diff=Math.abs(c1.radius-c2.radius)+Math.abs(c1.angle-c2.angle); 
				if (diff>1.0e-5 || c1.radius<0){
					System.out.println("********  l="+l+", m="+m+" phi:"+phi+ " ***********");
					System.out.println(c1+"=="+c2);
				}
				//System.out.println(c1); 
			}
		}
	}
	public static void test2(){
		float theta=2.3201294f;
		float phi=2.2299478f;
		test(theta,phi); 
	}
	public static void main(String[] args){
		test2(); 
	}
	public static void test(){
		
		for (int k=0; k<180;k++){
			float theta=k*PI/180.f; 
			for (int i=-180;i<180; i++){
				float phi=i*PI/180.f; 
				test(theta, phi);
			} 
		}
		
	}
}
